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We study isolated and binary neutron stars in dynamical Chern-Simons gravity. This theory 
modifies the Einstein-Hilbert action through the introduction of a dynamical scalar field coupled to 
the Pontryagin density. We here treat this theory as an effective model, working to leading order in 
the Chern-Simons coupling. We first construct isolated neutron star solutions in the slow-rotation 
expansion to quadratic order in spin. We find that isolated neutron stars acquire a scalar dipole 
charge that corrects its spin angular momentum to linear order in spin and corrects its mass and 
quadrupole moment to quadratic order in spin, as measured by an observer at spatial infinity. We 
then consider neutron stars binaries that are widely separated and solve for their orbital evolution 
in this modified theory. We find that the evolution of post-Keplerian parameters is modified, with 
the rate of periastron advance being the dominant correction at first post-Newtonian order. We 
conclude by applying these results to observed pulsars with the aim to place constraints on dynamical 
Chern-Simons gravity. We find that the modifications to the observed mass are degenerate with 
the neutron star equation of state, which prevents us from testing the theory with the inferred 
mass of the millisecond pulsar J1614-2230. We also find that the corrections to the post-Keplerian 
parameters are too small to be observable today even with data from the double binary pulsar 
J0737-3039. Our results suggest that pulsar observations are not currently capable of constraining 
dynamical Chern-Simons gravity, and thus, gravitational-wave observations may be the only path 
to a stringent constraint of this theory in the imminent future. 

PACS numbers: 04.30.-w,04.50.Kd,04.25.-g,97.60. Jd 
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I. INTRODUCTION 



Current astrophysical observations suggest that Gen- 
eral Relativity (GR) may have to be modified on large 
scales. Dark energy, dark matter and even inflation have 
been suggested as natural consequences of certain mod- 
ified gravity theories (see e.g. Ref. Perhaps similar 
modifications will be necessary in the non-linear, dynam- 
ical regime, the so-called strong-field, such as when black 
holes (BHs) merge and neutron stars (NSs) collide. In 
this regime, physical phenomena is not well-described by 
a leading-order weak-field and slow-motion expansion of 
the Einstein equations. Instead, one must either retain 
high-orders in such pcrturbative expansion, or solve the 
full set of Einstein equations numerically. The only way 
to determine whether modifications to GR in the strong 
field are necessary is to make observations and test GR 
in this regime. 

Modified gravity theories have been tested accurately 
in the Solar System and with binary pulsar observa- 
tions ■ While the former allow for tests in the weak- 
field regime only, the latter has led to tests when gravi- 
tational fields are moderately strong. Binary pulsar ob- 
servations, however, are still not capable of probing the 
non-linear regime of GR. For example, the orbital ve- 
locity of the double binary pulsar J0737-3039 [B-Q is 
roughly 10~^ smaller than the speed of light, and thus, 
its orbital behavior can be well-approximated by a weak- 
field, slow-velocity expansion of the Einstein equations. 



retaining only the leading and first subleading terms. On 
the other hand, GR will soon be tested during compact 
binary late coalescence through gravitational wave (GW) 
observations [l, [3j , which will allow for a probe of the 
fully non-linear and dynamical, strong-field regime. 

Recently, there has been a suggestion that binary pul- 
sar observations may constrain deviations from GR to 
such a degree that GW tests will not be necessary in 
the future 0. Such a suggestion emerged from studies 
of certain scalar-tensor theories [l^l, the constraints on 
which indeed cannot be improved with second-generation 
GW detectors, such as Advanced LIGO [uHli- But of 
course, this suggestion depends strongly on the particu- 
lar modified gravity model considered. For example, the 
same type of GW detectors may be able to place con- 
straints on massive graviton propagation that are three 
orders of magnitude stronger than current binary pulsar 
constraints |16l - [23 |. 

By studying NSs in dynamical Chern-Simons (CS) 
gravity [2^ , we here find another counterexample to this 
suggestion. Currently, the only constraints on this the- 
ory come from Solar System observations [l^l, by com- 
paring the CS correction to gravitomagnetic precession 
to observations with Gravity Probe B, and from table 
top [2^ experiments, by requiring that no CS correc- 
tions are present above the smallest gravitational length 
scales sampled experimentally on Earth. These tests lead 
to comparable constraints, namely -y/a < 10^ km, where 
a is the dimensional CS coupling constant. Such an in- 
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credibly weak constraint is perhaps not surprising, given 
the weakness of the gravitational field in the Solar Sys- 
tem. Recent work has shown that future GW observa- 
tions of BH binary inspirals could improve upon Solar 
System constraints by as much as seven orders of mag- 
nitude j26l - l29l |. In this paper, we study whether current 
isolated and binary pulsar observations can already con- 
strain dynamical CS gravity. We will see that indeed 
this is not possible and that only future GW observa- 
tions of compact binary coalescence can place stringent 
constraints on this theory. 

Dynamical CS modified gravity is a parity-violating, 
quadratic-curvature theory that is defined by modifying 
the Einstein-Hilbert action with a term that is the prod- 
uct of a dynamical scalar field and the Pontryagin den- 
sity (contraction of the Riemann tensor and its dual). 
The scalar field has dynamics through the addition of a 
standard kinetic term and a potential to the action. The 
Pontryagin term in the action is involved in anomaly can- 
cellation [soj . Such a term also appears naturally in het- 
erotic superstring theory [30l - l32| and loop quantum grav- 
ity jssl . [3^ . Dynamical CS gravity also arises naturally 
in effective field theories of inflation [3^1 • Historically, 
CS gravity was introduced without scalar field dynam- 
ics, assuming the field was given a priori [s^. Such a 
formulation was plagued with problems, such as metric 
instability [37| and overconstrained field equations. The 
modern (dynamical) incarnation of the theory restores 
the dynamics of the scalar field through a standard ki- 
netic term in the action (38| , while treating the model as 
an effective theory, thus avoiding instabilities (soj . 

Dynamical CS gravity has only recently begun to be 
studied in detail. Non-rotating BHs are described by the 
Schwarzschild metric, but rotating ones are not. Ana- 
lytic slowly-rotating BH solutions have been constructed 
to linear order in spin [s^, and to quadratic order 
in spin (25j within the small-coupling approximation, i.e. 
linearizing all expressions in the CS coupling. Slowly- 
rotating NS solutions were first constructed in [4l| to 
linear order in spin and within the small- coupling ap- 
proximation. This work was extended in |24l |. where NS 
solutions were obtained still to linear order in spin, but 
without imposing the small-coupling approximation. 

The purpose of this paper is to investigate whether 
any meaningful constraints on dynamical CS gravity can 
be obtained from isolated [l^ and binary pulsar obser- 
vations In order to achieve this goal, one must first 
study isolated NSs solutions to quadratic-order in spin, 
since the latter will introduce modifications to the NS 
mass as measured by an observer at spatial infinity. One 
must then place such NS solutions in binary systems and 
study the conservative and dissipative corrections to the 
evolution of the binary. The former come from the CS 
deformation of the quadrupole moment of each individual 
NS, as well as from the scalar dipole-dipole interaction 
of the binary components. The latter are caused by the 
modification to the emitted gravitational and scalar ra- 
diation. These corrections can be calculated once one 



finds the scalar dipole charge induced by each individ- 
ual NS and its associated CS quadrupole moment defor- 
mation. The charge can be extracted by studying the 
asymptotic behavior of the scalar field at spatial infinity 
to linear order in spin. The quadrupole moment deforma- 
tion appears at quadratic-order in spin. All throughout 
this paper, we treat dynamical CS gravity as an effective 
field theory. Such a treatment is definitely valid for the 
systems considered here, as explained in detail in [25| . 
In turn, this implies that we can use the weak-coupling 
approximation to linearize all expression in the CS cou- 
pling constant. Such a treatment guarantees that the 
field equations are second order, and the theory is ghost- 
free Hi. 



A. Executive Summary 

The first half of this paper focuses on finding iso- 
lated NS solutions to quadratic order in spin in the 
small-coupling approximation. We follow Hartle's ap- 
proach [43|, in which one treats rotation perturbatively, 
i.e. one expands in powers of the product of the NS mass 
and the NS rotational angular frequency. As expected, 
we find that CS corrections at quadratic order in spin 
modify the NS mass monopole and quadrupole moments, 
where the former can be absorbed by a redefinition of 
mass. At linear order in spin, the CS corrections ap- 
pear at much higher multipole order than quadrupole, 
since the correction to the current dipole moment can be 
absorbed by a redefinition of spin. Therefore, CS cor- 
rections at quadratic order in spin can be larger than 
those at linear order if NSs are spinning moderately fast, 
yet sufficiently slowly for the small-rotation expansion to 
hold. 

Figure [1] shows the scalar dipole charge p, (top left) and 
the CS corrections to the mass SM (bottom left), spin 
angular momentum SJ (bottom right) and quadrupole 
moment Q (top right) as a function of the GR mass pa- 
rameter in solar masses Mq. In this figure, C is the 
dimensionlcss coupling constant of dynamical CS grav- 
itjQ, ^* is the NS angular velocity, films is the angu- 
lar frequency of a NS with a period of 1ms, J is the 
NS spin angular momentum in GR and we used the 
APR H, SLvim^, Lattimer-Swesty (LS220) [13, 13 
and Shen j48l - l5G | equations of state (EoSs). Observe that 
the CS corrections reduce the observed mass, but they in- 
crease the observed angular momentum and quadrupole 
moment with increasing M*. 

Such corrected NS observables can then be contrasted 
with NS observations to constrain dynamical CS gravity. 
Observations of the massive millisecond pulsar J1614- 
2230 [42} require that one be able to construct NSs with 



^ This quantity is proportional to the fourth power of the Chern- 
Simons natural length scale and it is such that as f — > one 
recovers GR. See Eq. l(9]| for a more precise definition. 
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observed masses larger than 1.93Mq. Unfortunately, but 
perhaps not surprisingly, the allowed maximum NS mass 
in dynamical CS gravity depends not only on the EoS but 
also on the magnitude of the dimensionless CS coupling 
parameter ( oc ^csM^ /T^t- Figure [2] shows the mass- 
radius relation for different EoSs with different choices of 

Therefore, such a test of dynamical CS gravity cannot 
be performed due to degeneracies with the EoS. 

The CS corrections to the NS quadrupolc moment and 
the moment of inertia induce modifications to the or- 
bital evolution of binary pulsars, which then can be con- 
trasted with binary pulsar observations to constrain the 
modified theory. Observations of the double binary pul- 
sar J0737-3039 0-0] require that the evolution of certain 
Keplerian elements 0, 0, HH agree with GR to within a 
certain observational uncertainty. The CS modification 
to the quadrupole moment leads to corrections to these 
evolution equations, which we explicitly calculate in this 
paper. We find that the rate of periastron advance w is 
the best observable to constrain CS gravity, because its 
CS modification enters at first PN order relative to the 
GR behavior, i.e. this modification is suppressed by a 
single factor of 0{m/a), where m is the total mass of 
the binary and a is the semi-major axis of the binary. 
However, for J0737-3039, this still implies a suppression 
of the CS correction of 0(10"*) relative to the measured 
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FIG. 1. (Color online) Numerical (symbols) and fitted (curve) 
results as functions of NS mass with various EoSs. We plot 
the dimensionless scalar dipole susceptibility p, (top left), the 
quadrupole correction Q (top right), the mass shift SM (bot- 
tom left), and the angular momentum shift 5 J (bottom right), 
which are defined in Eqs. (|69|) . H66|) . and (jSSp respectively. 
The y-axis in the top-right panel is normalized by J'^/M, , 
the absolute value of the NS quadrupole moment in the point 
particle limit in GR, while the y-axis of the bottom-left panel 
is normalized by films, the angular frequency of a NS with 
a period of 1ms defined by (f2ims = 27r/lms). Moreover, Q, 
SM and 5J axe also re-scaled to = 0.1. Observe how the 
mass shift is always negative, while the quadrupole moment 
deformation is always positive. 
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FIG. 2. (Color online) Mass-radius relations for GR (solid 
curves) and CS with different coupling strengths {( = 
0.1,0.5,1 curves are respectively dotted, dashed, and dot- 
dashed). Each panel uses a different EoS: APR (top left), 
SLy (top right), LS220 (bottom left) and Shen (bottom right). 
The horizontal dashed line corresponds to the lower bound 
on the NS mass, 1.93M0, provided by observations of PSR 
J1614-2230. The mass-radius relations depend on the spin of 
the NS, which we have set to match that for PSR J1614-2230. 
Since APR, SLy and Shen EoSs can each produce NSs above 
the mass bound for all < 1, there is a clear degeneracy be- 
tween EoS and the CS correction to the mass-radius relation. 
Thus, one cannot constrain the theory from observations of 
PSR J1614-2230. 



GR Lo. Moreover, this CS effect is even smaller than the 
GR spin-orbit correction to w, which enters at 0.5PN or- 
der and depends on the NS's moments of inertia. There- 
fore, current binary pulsar observations arc not accurate 
enough to allow for tests of dynamical CS gravity. 

Unlike for scalar-tensor theories [1], our results sug- 
gest that dynamical CS gravity cannot be constrained 
well with binary pulsar observations using current data. 
Instead; GW observations may be the only way to place 
stringent constraints. We have therefore found a con- 
crete example of a modified gravity theory that cannot 
be stringently constrained with current electromagnetic 
observations. 



B. Organization 

The organization of the paper is as follows. In Sec. |lll 
we introduce the basics of dynamical CS gravity and ex- 
plain the approximations that we use throughout the pa- 
per. In Sec. mil we explain the coordinate systems that 
we use and impose an ansatz on the metric and the stress- 
energy momentum tensor of the matter field. In Sec. IIV| 
we derive both GR and CS field equations at each or- 
der in spin and derive the exterior solutions, modulo in- 
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tegration constants, which are discussed in Sec. [V] In 
Sec. I VII we explain how to construct the interior solu- 
tion. This amounts to finding the matching conditions 
at the NS surface and the initial conditions at the center. 
Then, we explain the numerical procedure that we use 
to solve these sets of differential equations. In Sec. IVIH 
we present the numerical results obtained by solving the 
interior equations and matching the solutions to the ex- 
terior ones. We derive the evolution of the NS binary 
in Sec. I Villi and in Sec. |IX] we apply the results to the 
massive millisecond pulsar J1614-2230 and to the double 
binary pulsar PSR J0737-3039. We conclude in Sec. H 
with a discussion of possible avenues for future work. 

All throughout the paper, we mostly follow the con- 
ventions of Misner, Thorne and Wheeler [s^]- We use 
the Greek letters (a, /3, • • • ) to denote spacetime in- 
dices. The metric is denoted by g^ti^ and it has signature 
(— , +, +, +). We use geometric units, with G = c = 1. 



II. DYNAMICAL CHERN-SIMONS GRAVITY 
A. Basics 

The action in dynamical CS gravity is given by (23j 

a 



By using this equation, one can show that 




(1) 



Here, Kg = (167r)^^, g denotes the determinant of the 
metric and Rf_i,^pa is the Riemann tensor. *Jlf^'^P'^ is 
the dual of the Riemann tensor which is defined by [2^ 



(2) 



where e^"^"^ is the Levi-Civita tensor, is a scalar field 
while a and (3 are coupling constants. £„,at denotes the 
matter Lagrangian density. Following Refs. [25l. Isol l53j. 
we have neglected the potential of the scalar field for 
simplicity. We take and /? to be dimensionless while 
we set a to have a dimension of (length)^ [39j . 

The field equations in dynamical CS gravity are given 
by HI 



ur ^ Our ^ A*^ ^fiu) ' 



(3) 



where G^i/ is the Einstein tensor and Tj^J^^ is the matter 
stress-energy tensor. The C-tensor and the stress-energy 
tensor for the scalar field are defined by 

C'' = iy o.R'\ + {V s^T R^^^"^" , (4) 

(5) 



pu 



The evolution equation of the scalar field is given by 

a 



4/3 



(6) 



1 



(7) 



Together with the Bianchi identity, if we take the diver- 
gence of Eq. ([3]) , we end up with 



(8) 



Thus, the equation of motion for a test particle is not 
modified. 

The evolution equation for the scalar field, Eq. (jS]), ad- 
mits a flat background solution of the form d = C^a;^, 
for ~ const, i.e. a solution to the homogeneous equa- 
tion in a fiat background. The constants control the 
strength of the CS modification, and thus, one would 
expect them to be proportional to a/ j3. Although such 
a solution is in principle allowed, would have to be 
prescribed a priori^ as some form of cosmological term. 
Moreover, such a term would ruins the 'd shift-invariance 
of the equations of motion, which is a desirable prop- 
erty that allows us to treat d as massless. Also, such a 
scalar field would have an infinite energy density, when 
its stress-energy tensor is integrated over the entire man- 
ifold. Finally, dynamical CS gravity with such a homo- 
geneous scalar field is very similar to the non-dynamical 
theory, and thus, it is already severely constrained by cos- 
mological observations 5J]. For the rest of this paper, we 
will ignore the homogeneous solution for -d, and instead, 
concentrate on dynamically generated scalar fields. 



B. Small-Coupling, Slow- Rotation Approximations 

We here work in the small-coupling and slow-rotation 
approximations. The former implies that we treat the 
action given in Eq. ([1]) as defining an effective theory, 
which requires the CS quadratic term (the second term) 
to be much smaller than the Einstein-Hilbert term. For 
isolated NSs, the small-coupling approximation is valid 
if the following inequality holds: 



c 



7^6 



a 

(3Kg 



(9) 



where and 7?.* are the mass and radius of the NS 
respectively. ^^^^ correspond s to the characteristic length 
scale of the theory while ^JTZ^/M^ corresponds to the 
curvature length scale of the NS. This definition of the 
dimensionless small coupling constant C, is the same as 
defined in [23| modulo an 0(1) numerical factor. Figure 7 
of [23| shows that the small coupling approximation is 
valid when C 1. Solar System experiments require 
Ccs^ ^ C'(10^)km [23]. In this paper, we work to linear 
order in C^. 

Dynamical CS gravity must thus have cut-off scale be- 
yond which one cannot treat it as an effective model. 
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Ref. [25| estimated this scale by calculating the critical 
length scale at which loop corrections to the second term 
in Eq. ([T|) due to n-point interactions become of order 
unity. Requiring that this length scale be smaller than 
the smallest scales probed by table-top experiments led 
to the requirement that ^/a. < 0(10^ km). This require- 
ment is of the same order as the current constraints from 
Solar System experiments. For the NS systems we will 
consider here, both this requirement and the small cou- 
pling condition C ^ 1 are satisfied. 

We further assume that NSs are slowly rotating (i.e. 
Ixl ^ 1 where x = J/M"^ is the dimensionless spin pa- 
rameter, with J the magnitude of the NS spin angular 
momentum and M is its mass) and consider terms up to 
quadratic order. 

All physical quantities. A, can be expanded bivariately 



A 



(m,n) 



(10) 



where x' and a' are book keeping parameters to label the 
order of the slow- rotation and the small- coupling approx- 
imations, respectively. Notice that A(m,n) x^a". For 
the spherically symmetric case, Ru^ipcr* R'^'^'"^ = 0, and 
thus, there is no source to the inhomogeneous equation 
of motion for d. Therefore, for the scalar field, i?(o,n) ~ 
and 



^ = a' x'^ (1,1) +0{a'x'n 



(11) 



Notice that § 



(2m, 1) 



due to parity. 



III. SPACETIME AND MATTER 
DECOMPOSITION 

A. Metric 

Following Hartle [1^, we start with the metric ansatz 
given by 



s''('') {l + 2h{r,e))dt^ 
Mr) (i ^ ^rn{r,e) ^ 



r - 2M (r) 

+ {1 + 2k{r, 6)) IdO^ + sin^ 9 {d<j) ~ w(r, e)dt) 



(12) 



where D and A are 0{a"^x"^) quantities that only depend 
on r while h, k, fh and ui refer to perturbations that 
depencQ on both r and 9. M is related to A by 



M{r) 



1 _ e-^f'-; 



(13) 



These quantities do not depend on t and </> because we are search- 
ing for stationary and axisymmetric NS solutions. 



The coordinates {t, r, 9, (/>) are Hartle-Thorne coordi- 
nates. In particular, (r, 0) denote ordinary polar coor- 
dinates. 

Since we are treating rotation perturbatively, one must 
be careful with the choice of polar coordinates [1^. A 
perturbative approach is valid only when the fractional 
change in a quantity between rotating and non-rotating 
cases is small. If one were to use (r, 9) coordinates, this 
condition could not be met for the density p and the 
pressure p near the surface. This is because the rotation 
changes the shape of a star, and hence there are points 
on the NS surface where these quantities vanish in the 
non-rotating case while they acquire finite values in the 
rotating case, leading to infinite fractional changes. In- 
stead, we adopt the coordinates {R, 8) as proposed by 
Hartle l4j, which are defined via 



p[r(i?,e),e] = p(i?) = p(o.o)(i?), &^o. 



(14) 



In other words, the new radial coordinate R is chosen 
such that the density at p[r{R,Q),Q] in the rotating 
configuration is the same as p(o^Q){R) in the non-rotating 
configuration. By construction, the density and the pres- 
sure in the new coordinates contain the non-rotating part 
only: 

p{R) = P(o,o){R), piR) ^ P(oM^) ■ (15) 
We expand r{R, 9) as 

r{R,Q) = R + ^{R,e), (16) 

where 

^R, 9) = a'2x"e(2,2) {R, e) + 0{a'\'^a'\'^) . (17) 

Notice that we have neglected 0{a"^x'^) quantities (pure 
GR, quadratic in spin effects) in Eq. p7)) . Since in this 
work we are interested in the detectability of the CS cor- 
rections, and the mentioned terms would simply add lin- 
early but not modify the CS corrections at ©(a'^x'^): 
ignore the 0{a"^x'^) term henceforth. 

After the coordinate transformation, the new metric is 
given by 



l + 2h + ^ 



. dv 
dR 



- R^uj^ sin^ 9 



dt^ 



-2R^Lo sir? <ddtd(t) + [i?2(l + 2k) + 2R(\ sit? Qdcj)^ 



2m 



R - 2M 



_dX 
dR 



+2e^^dRde + [i?2(l + 2k) + 2R^] dQ^ 
a9 

+O(a'0x'^a'\''), 



where 



M{R) 



(l-e-^(«))i? 



(18) 



(19) 
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Each quantity in the above equation is bivariately ex- 
panded as 



A(i?) 

uj{R,e) 

h{R,e) 
m{R,e) 
k{R, 6) 



:x'^(i,o)(-R,e) + a'2x^i,2)(i?,e) 
^a''x''^2,2)(i^,e) + o(«'°x'^a'V'), 

-a'\''m^,.2){R.Q) + 0{a"x'^a'\'')., 

^«''x"fc(2,2)(i?,e) + o(«'V^«'V). 



(20) 



Due to parity, the (i, t), (i?, R), (8, 9), (0, ^) and (R, Q) 
metric components only have terms proportional to even 
powers of %' while (<, 0) component only contains terms 
that arc odd powers in x' ■ Also, 0{a') terms do not 
appear in the metric because of the structure of the field 
equation and the fact that is proportional to a in the 
small-coupling approximation. 

Notice that the coordinate deformation ^{R, 8) is well 
defined only inside the star. We take it to be constant 
outside the star. This means that the exterior metric in 
(i, r, 0, (j)) coordinates can be obtained simply by replac- 
ing R r and — > 6* in the exterior metric in {t, R, 0, 4>) 
coordinates. 



B. Neutron Star Stress-Energy Tensor 

In this paper, we assume that the matter field inside 
a NS is a perfect fluid and that the NS is rotating uni- 
formly. The stress-energy tensor for the matter field T^'"' 
is given by 



= {p + v)u^,Uy + pg^^ , 
where the four-velocity is given by 
u'' = (u°,0,0,rj,w°) . 



(21) 



(22) 



is the constant angular velocity of the NS. By using 
the normalization condition u^u^ = —1, we obtain the 
time component of the four-velocity as 



dv 



1(2,2) 



(23) 



As mentioned in Eq. (|5]), the stress-energy tensor of 
this minimally-coupled fluid is divergence free, which 
contributes another equation to the system. In order to 
close the system, one needs a relationship between den- 
sity p and pressure P, an EoS. In this paper we only 
consider one-parameter EoSs, i.e. P = P{p) without any 
entropy dependence. The EoSs we use are APR [44l |. 
SLy |45l . kg , Lattimer-Swesty with nuclear inc omp ress- 
ibility of 220MeV (LS220) [43, SI and Shen For 



the latter two, we use a temperature of O.OlMeV and an 
electron fraction of 30%. The EoS and conservation of 
matter stress-energy tensor close the system. We discuss 
the matter equations of motion order-by-order below. 



IV. MODIFIED FIELD EQUATIONS 

In this section, we derive the modified field equations 
expanded in x' a-nd a' . 



A. Zeroth Order in Spin 

First, we consider equations at 0(x'°). As mentioned 
previously, in the spherically symmetric case there is no 
CS correction, and hence the field equations reduce to 
the Einstein equations for a non-rotating star. The only 
non-vanishing components are the (t,t), {R,R), (0,0) 
and (0, (j)) ones, but the last two components are linearly 
dependent. The first two components yield 



dM , 
dR 



^ AttR^P + M 
^ R{R-2M) 



(24) 
(25) 



Together with the EoS, we need one more equation to 
close the system of differential equations. Instead of using 
the (0, 0) or ((/i, (j)) component of the Einstein equations, 
one can use the R component of the conservation equa- 
tions of the matter stress-energy tensor [Eq. (|5])] which 
is given by 



dv 
dR 



dp 



p + p dR 



(26) 



By combining Eqs. (P5|) and ([2^. we obtain the Tolman- 
Oppenheimer-Volkoff (TOV) equation 



(27) 



dp _ {4:TtR^p + M){p + p) 
dR ^ RiR - 2M) 

B. First Order in Spin 



1. GR 

At ©(a'^xOi the only non- vanishing component of the 
Einstein equations is the {t, (p) one. By using Eqs. (p4)) - 
(P7|. the associated equation can be simplified to 



(1.0) 



1 - iTR^{p + p)e^ 5w(i,o) 



902 

167r(p 4- p)e'*'a)(i_ 



R 
3 cot© 



dR 



(1,0) 



de 



0) 



(28) 
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where 

^(1,0) = f^* - ^^(1,0) • (29) 

We simplify the equation further by performing a Legen- 
dre decomposition: 

d 



1 



Pc(cose) , (30) 



sineae 



where Pe is the ^-th Legendre polynomial. Then, Eq. ([28 
becomes 



di?2 R 
-(^ + 2)(£-l) , 



i?2 



dR 



0. (31) 



By imposing asymptotic flatness and regularity at the 
center, one can show that uj^ must vanish for all £ except 
£ = 1 [43] • Thus, the above equation reduces to 

d^wi I- 'KR^{p+p)e^ duJi 



- 167r(p + p)p/aji = 



R 



dR 



(32) 



2. CS: Scalar Evolution Equation 



Next, we look at the scalar evolution equation at 
0{a'x')- By using Eqs. ([Ml) -(HZ]), the equation can be 
written as 



cote 



R 

5^(1,1) 



9i? 



99 



d 



2,~, 



8^1^5 e(A-.)/2f 31119^1;^ + 2 cos e^'^^^'") 



where 



dRBG 



M 



dR 



(33) 



(34) 



(4/3)7ri?3 

denotes the density shift from the average. Since we are 
only interested in stationary and axisymmetric solutions, 
we consider a "d that only depends on the coordinates R 
and 6. As in Eq. (pO| . we perform a Legendre decompo- 
sition of -d: 



i9f(i?)P£(cose), 



(35) 



and find that the £ = 1 term of Eq. ()33p is the only one 
with a source term. The only homogeneous solution that 
satisfies both asymptotic flatness and regularity at the 
center is the trivial i?^ = for £ 7^ 1. For £ = 1, Eq. dSS]) 
becomes 



i? dR R'^ ^ 



dR^ 

B dR 



(36) 



5. CS: Field Equation 

Next, we consider the modified Einstein equation at 
0{a'^x')- in GR, the only non- vanishing component 
of the field equations is the {t, (j)) one. Notice that since 
■i? is of C(a'x'), there is no 0{a'^x') contribution from 
T^^. The (t, (j)) component of the field equations is given 
by 



9i?2 



W(i,2) 



e 

P2 \^ 992 

128^2^g(.+A)/2 

i?3 sin e 



5i? 



i? 

3 cote 
92^(1 



OR 



dujn 



(1,2) 



9e 



1,2) 



9i?9e 



5^(1,1 



9e 



As for W(i,o)j 'we decompose W(i,2) via- 

1 d 



(37) 



^(i,2)(i?,e)-^u;,(i?) 



sinese 



Pf(cose) . (38) 



Again, we find that Eq. ([57)1 with £ 7^ 1 becomes a ho- 
mogeneous equation, leading to tu^ = for £ 7^ 1 after 
imposing asymptotic flatness and regularity at the cen- 
ter. For £ = 1, the equation reduces to 

d'^wi l-TTR^{p+p)e^ dwi ^ 

lbTr[p + p)e wi 



dR"^ R 

128^2^g(.+A)/2 



dR 



i?3 



S R 



d^i 
dR 



^1 



(39) 

Notice that Eq. (f32|) corresponds to the homogeneous 
version of Eq. ([39]) . 



C. Second Order in Spin 

In this subsection, wc derive the stress-energy conser- 
vation and the fleld equations at quadratic order in spin. 
As mentioned previously, we do not consider the GR part 
since 0{a'°x'^) terms do not appear in any equations at 
0{a'^x'^)- Following the previous subsection, we per- 
form the Legendre decompositions: 

/i(2,2)(i?,e) = Ef/i£(^)^f(cose), (40a) 
m(2,2)(i?,e) = m^(i?)Pf (cose), (40b) 
fc(2,2)(i?,e) = j:fkiiR)Peicose), (40c) 

C(2,2)(i?,e) = ZiUR)Piicose) . (40d) 

We use the gauge freedom of the theory to set fco(P) = 

M- 



1. CS: Stress-Energy Conservation 

The non-vanishing components of the stress-energy 
conservation equations V^T^'J"' = are v — R and e. 
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First, let us look at the component. By performing a 
Legendre decomposition, one finds that ^ = 2 is the only 
non-vanishing mode. By using the previously obtained 
results, one arrives at the algebraic condition 



6 



R{R-2M) 
3(47ri?3p + M) 



(3/i2 - ie-^R^QiWi) (41) 



in the interior of the star. 

Next, let us look at the R component. The only non- 
vanishing modes are the ^ = and H = 2 ones. For the 
latter, we obtain Eq. (PT|) after we integrate the differ- 
ential equation once with respect to R. For the former, 
again by integrating once, we obtain the algebraic rela- 
tion 



AnR^p + M 2 



R^ujiWi + hoc 



(42) 



valid in the interior of the star, where /iqc is an integration 
constant that corresponds to h{R = 0). 



2. CS: Field Equations 

Let us now derive the modified field equations at 
©(a'^x'^)- The only non- vanishing modes are the i = 
and £ = 2 ones. Let us first focus on the £ = case, where 
our dependent variables are /iq, mo and ^o- Eq. (j42|) al- 
ready determines /iq, and thus, we are left with 2 degrees 
of freedom (mo and ^o)- Following Ref. [i^, we use the 
(t, t) and (i?, R) components of the field equations, and 
Eq. (|42l), to find 



dniQ 
~dR 



dR 



6 



dwi \ dull 



^-^e-^RHp + p)^iWi+''=^l3er'R' 



dR 
ddi 
'dR 



dR J dR 



An 



An 



^g-(.+A)/2 ) ^ (g^-A ^ ^ ^ 8^^2p) _ 2 [2e-^ + 3 + 16nR\p + p)] 



dR dR 



dwi 
'dR 



647ri?2 (p+p)cZ,^^_327ri? 
ait 



, dp 



^J^ + iP + P) [2-{l + SnR^p)e^] 



(43) 



d^o _ 1 f 6(1 + SnR^p) 

'dR ~ 3[1- (l-|-87ri?2p)e^] \ 'r 



e^^mo -I- 3- 



X^nR^p - (1 - 87ri?2p)(i + %T:R^p)t 

'r 



dR dR 
2[8^i?2(p + p)-l]^i^ 



^ dwi 



4 (3e"^ - 1 - ?>nR^p) R^e-'^'^'^^CjiWi + AnpR^ 



ddi 
'dR 



(3, 



,-A 



1 - SnR^p) R 



d-d I duj 



dR dR 



\ + l6ne^{p + p)diCoi 



(44) 



Using Eqs. (021), (|43]) and (04]), one can also obtain 



dho 

Ir 



SnR^p + 1 



-e^^mo 



4nR'^{p + p) 1 3 8ne-''R'^{p+p) ^ 



-ho + -R'^e' 

R-2M 6 dR dR 3 R{R-2M) 



An „e 



2n 



3 ' R 



ddi 



3e-^ - 1 - SnR'p) ( + 167re^p + p)^^^;^ ] + 2 



^8n R^{p + p)^l d^i 
R ^ dR 



dR, 
(45) 



We now focus on the i ~ 2 mode, where we have 4 unknown functions: /12, m2, k2 and ^2- Eq. (|4ip already gives 
us ^2, and thus, we need 3 more equations to close the system. We choose one of these to be the (0, 9) component 
of the field equations minus the (0, 0) component. This yields 



m2 



-e~^Rho 



-Pe-^Rdl 



g-(,.+A)^5 



dwi duj 



dR dR 



i + 16n{p + p)wiu}i 



16n 



^g-(.+3A)/2^ 



e~^R 



duJi 



ddi 



[e-^ + AnR^ [p + p)] - STri?^ (p + p)Cj^ _i _ ( i?^ - {AnR^p + M) {p + p)] diCji 



dR 



dp 



dR 



d-di dijJi 
'dR~dR 

(46) 
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For the remaining 2 equations, wc use the {R,Q) and {R, R) components of the field equations. By using Eq. (|4ip . 
they yield 



dk2 



dh2 (3e-^ - 1 - SnR'^p) 



dR 



2R 



(47) 



and 



[3-47ri?2(p + p)]e> 



1 + (1 + 8TTR^p)e^ dk2 
2 



2 — K2 H 



3 ^ i? 



"6~ dR dR 



in 



i?2 

.(A-.)/2 ^3g-A 



7712 



47r ZdT?! 



1 - SttR^^p 



~3~ 



2 ^ d7?i dwi 



it ait 3 



(48) 



These equations are technically valid only inside the NS, because they were developed from the matter conservation 
equations [Eqs. (PT|) - (|^^ ]. However, we have checked that the field equations outside the star are identical to those 
above after setting p — — p. 

I 



ISOLATED NEUTRON STAR SOLUTION: 
EXTERIOR FIELDS 



In this section, we analytically solve the equations 
derived in the previous section outside the NS, where 
p = = p. Wc use the superscript "cxt" to refer to exte- 
rior quantities. We impose asymptotic flatness at spatial 
infinity 0jS cL boundary condition. 



A. Exterior Solutions 

1. GR 

At 0(a'°x'°), we solve Eqs. ^ and ^ to obtain 

2Ai* 

(49) 



R 



where iVi* = const, is the mass of the non- rotating NS, 
which is to be determined by matching v^^'^ or A'^'^* with 
j,int ^int ^j^g surface. At 0(a'°x'^), we substitute 
the above equation in Eq. ([32]) and obtain 



2J 

i?3 ' 



(50) 



2. C'S: First Order in Spin 



Let us substitute Eqs. p9|) and (|50|) in the equation 
for 7?i [Eq. ([55)) ]. The solution to this equation is 

t 5a J 1 f_ „ili* . 18Af2 

^1 



8 P Ai2 i?2 

i?2 - 



1 



1 



R 
2MI 



R 



1 



5 i?2 



R 



\nf{R) 



(51) 



where is an integration constant that is to be deter- 
mined by matching this solution with the interior one 
at the NS surface under continuous and differentiablc 
boundary conditions. Notice that an integration constant 
has been set to by asymptotic flatness. 

The modifled field equation at ©(a'^x') can be solved 
by substituting -df^^ in Eq. ([39]) to obtain 

27 Mi 



2Ja 



5 CcsJ 



R^ 8 il/*i?6 



R 



+ 



R 
2MI 



12 AL 

54 1/3 

48 Mf ■ 
"5""r4 



10 i?2 



64 



(52) 



Here, Jcs is an integration constant that corresponds to 
a CS correction to the spin angular momentum (there 
is also a correction to the spin angular momentum from 
the log term). Eqs. ([5T]) and ([5^ agree with those found 

in m. 



where J is an integration constant that corresponds to 
the spin angular momentum of the NS in GR. We define 
the mass and the spin angular momentum by expanding 
the metric about R = oo and extracting the appropriate 
coefficients [24i. 



3. CS: Second Order in Spin 

Next, we find the exterior solutions at 0{a'^x''^). For 
the £ = mode, we can solve the exterior version of 
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Eqs. (031) and (gSj) to find 



768M3i?5/(i?) 
AMI 2 M 3 



1 



A/, 
24 ^ 



A/f 684 Af 3 m4 ^4 



1 



' R 



R 



5 i?4 2Af* 



A4 lA/2 IA43 34Af4 24A45,, ^, , 



15^2 
" 2 ^^A/4 

768 A/4i?3 
'Af2 

i?4 



1 



6A4 

.A//* 



1 



^)ln/(i?) 



■ 2C,9 



1 



322 A/2 
Ah 4 A/2 72 m3 



12A/2 
198 A/3 6276 A/f 



/(/?)[ln/(/?)]^ 

17496 A/f 



A/n 



2JJc 



/?3 

R 



A/f 



/? 

aT 



1 - 2 



Ai, 
/? 



2 /?2 



2A/. 
ln/(/?) 



175 /?4 

1 A^ 
5"r2" 

/?2 



1 + 



4A/2 



25 /?5 
78 A/3 

A/* 



A/'* , 
24-^)ln/(/?) 



2JJc 



+ A/c 



(53) 



(54) 



where A/cs is an integration constant that corresponds 
to a part of the correction to the mass if one expands 
the above solutions about R = 00. For the £ = 2 
mode, we obtain similar solutions, shown in Appendix 
but with another integration constant Cq. If one is to 
deal with a BH, one needs to impose regularity at the 
horizon. By setting Jcs — 0, C^} = 0, A/cs = and 
Cq = (3015/14336)Ccs(J^VA4^), the exterior metric in 
(t, r, 9, (p) coordinates reduces to the previously found BH 
solution 1251. 



B. Asymptotic Behavior at Spatial Infinity: 
Scalar Dipole Charge and mass 
Quadrupole Moment Deformation 



The asymptotic behaviors of hfj 
tial infinity is 



cxt 



and wf^^ about spa- 



7 ext 

"0 



with 



6M = A/cs + 



25J 

1^' 



25 Cp^ 

7SCS " 



1536 



A/J 



SJ^Jn 



25 Ci)J 
^^cs-^ 



(55) 

(56) 
(57) 



Physically, 6M and SJ correspond to the CS corrections 
in the mass and the angular momentum respectively. 

Let us now define the observable mass A/ and angu- 
lar momentum J as measured by an observer at spatial 
infinity via 



such that 



ext / -1 

9tt 1 - 



2A4 
R 



2J 



■ sm' 



•'9, (59) 



near spatial infinity. The exterior metric and the scalar 
field obtained in the previous subsection can be rewritten 
in terms of A/ and J by replacing A/* — > A/ — 6M and 
J-5J. 

The full exterior metric can then be written as g'^ = 
^cxt,GR_|_^oxt,cs^ where g'^'^^ is the Hartle-Thorne metric 

with mass M and angular momentum J, while g'j^''^^ is 
a CS correction. The latter, in (<, r, 9, (p) coordinates, has 
the following asymptotic behavior about spatial infinity 



ext.cs 

9tt 



1 25^csC'i3(C^ 


-8)J2h 


- 2048C<; 


,A/8 


3840 


A/5r3 






X (3 cos^ 61 - 1) + 


°(^) 






1 25^asC4C^, 


-8)J2h 


~ 2048Q 


,A/8 


3840 


A/5r3 






X (3cos2 6i- 1) + 








1 25^csC'i9(C'i9 


- 8)P H 


- 2048C(; 


,A/8 



3840 



X (3cos^6i- 1) + © 



A/5r 
1 



sin^ 9g° 

^Ccs 



(3-C^)J 
A/r4 



sm 



4 



(60) 

(61) 

(62) 
(63) 
(64) 



A/ EE A/* + 6M, J = J + 5J, 



(58) 



One can read off the correction to the gravitational po- 
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tential per unit mass from Eq. ((60|) as 



(65) 



with 



^_ 1 25^csC^iC^~S)P + 20A8CQM^ ^^^^ 



7680 



M5 



Here, 5* is the unit spin angular momentum vector of the 
NS and n* is the unit vector to the field point. Notice 
that this definition of S" differs from that used in [53j . 
because in the latter 5* was not a proper unit vector. 
Q corresponds to the correction to the quadrupole mo- 
menllfl and in the non-spinning BH case, Eq. (|66p re- 
duces to Qbh = (201/3584) (^cs/M'^)M3x^ which agrees 
with [H. 

Finally, the asymptotic behavior of the scalar field is 



5(3 — C^) a ,cos( 



24 



o 



1 



(67) 



where x = J/M"^. Comparing this to •§ = [i^nijr^ 
0(r~^^ (see Eq. (57) of Ref. [53) where /i* is the scalar 
dipole charge, we can extract /z': 



5(3 ~ Cjf) g 
24 



(68) 



For later convenience, we introduce the dimensionless 
scalar charge \i in terms of Cf) and the compactness of 
the NS C: 



3-Ctf 1 



C3 ' 7^, 



By using this /i, Eq. (|68|) can be rewritten as 



(69) 



(70) 



From Eq. (59) of Ref. [ll], we expect pJ- to be propor- 
tional to C'^, and hence, we have factored out from 
the definition of fi above. In the non-spinning BH case, 
fi reduces to /Ibh = 8, because then = 1/8, so that 
MbhC^ = 1. 

Notice that the difference between M and Af* (and 
also between J and J) in Eqs. (|60p - ([7D|) is irrelevant to 
the order of approximation considered here. Hence, we 
can freely set M = (and J = J) in these equations. 



^ This definition of quadrupole moment is difi^erent from the usual 
one given in e.g. [ssl . ISfiH by a factor of 2. Although somewhat 
unconventional, we choose here to continue using the definitions 
of Ref. fH. 



VI. ISOLATED NEUTRON STAR SOLUTION: 
INTERIOR FIELDS 

In order to determine /I, Q, bM and (5J, we need to 
calculate the integration constants that appear in the 
definitions of these quantities. This can be achieved by 
matching the interior solutions to the exterior ones at the 
NS surface. In this section, we first explain the matching 
conditions at the surface. We then obtain the asymptotic 
behavior of the metric perturbations at the NS center to 
obtain the initial conditions for the numerical integration 
of the interior solution. We conclude with an explanation 
of the numerical procedure adopted in this paper. We use 
the superscript "int" to refer to interior quantities. 



A. Boundary Conditions at the Surface 

In the previous subsection, we imposed asymptotic 
flatness at spatial infinity to find the exterior solutions 
but 4 integration constants remain, which are to be de- 
termined by the boundary conditions at the NS surface. 
At 0(a'"x'"), we solve Eqs. ([25]) and together 
with the EoS, for A, p and p. At the NS surface, we 
impose the continuity condition 



P 



2M* 



(71) 



The first equation determines TZ^ while the second equa- 
tion determines M*. 

At 0{a'°x'), we solve Eq. ([32]) for uji. Since this is 
a second-order differential equation, we need to impose 
two boundary conditions at the surface. We choose one 
to be the continuity condition for i.e. 



2 J 



(72) 



The other condition can be obtained by integrating 
Eq. (j32p from 7^* — e to 7?,* + e and taking the limit 



as e 



dR 





r dcjf^ 


= lim 




e-i-O 


[ dR 






= lim 












~ lim 









(TZ, + e) 



dR 



(7^, 



d'^uji 



dR, 



1 - TiR^{p + p)e^ dull 



R 

+167r(p + p)e^wi} di?, 
where we have introduced the notation 
[A{n,)] = lim [^°'='(7^* + e) - A"'\n, 



dR 
(73) 

0] (74) 



for any quantity A(R). Since the integrand is bounded, 
the right hand side of Eq. ()73|) vanishes, leading to 



-dR^""* 



= 0. 



(75) 
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Let us introduce the moment of inertia /, defined 
by El 



Equation ((72|) can then be rewritten as 

21 



(76) 



(77) 



By Eqs. (HH), ((25|. (f32| and ([75]), / can be expressed in 
integral form as [43, |57| 



3 r!, i?-2A/ 



(78) 



This expression is valid both in GR, as well as in dynam- 
ical CS gravity H. 

As in the ©(a'^x') case, the scalar evolution equation 
[Eq. (p6l) ] is a second-order differential equation. Since 
all terms in this equation are bounded at the surface, we 
impose the continuity and differentiability conditions 



[i?i(7e.)] = o = 



(79) 



At ©(a'^xOi the non- vanishing component of the field 
equations is given by Eq. pop . which is also a second- 
order differential equation, but here the situation is 
slightly different. First, we impose the continuity con- 
dition at the surface 



Next, we integrate Eq. 
take the limit e — >■ as 



K(7^*)] = 0. 

from 72.* 



e to 72.* 



(80) 
e and 



lini 



lim 



dR, 



R 



-167r(p -I- p)e^wi 
1287r2ae('-+^)/2 

i?3 



1 - nR'^ip + p)e^ dwi 
dp 



SR 



ddi 
liR 



R 



dR 



dR, 



on,+e g(.+ A)/2^ ^ 

-1287r2alim / -^dR. 



n,-t 



i?2 dR 



where in the last equality, we only kept the term in the 
integrand that is not boundecQ. This equation can be 



further simplified by integrating by parts: 



= -12877^0; lim 



e->0 



TZ,-e 



dR 



^2 



dR 
dR 



^2 ' 



(82) 



Since the integrand in the second term is bounded, this 
term vanishes, so that one finally finds the jump condi- 
tion 



^(72*) 
dR ^ *' 



U8n'a^^p{n.). (83) 



This agrees with the condition found in 12J] in the limit 
e — >■ 0. Since the density of the NS at the surface is typ- 
ically about 7 orders of magnitude lower than the mean 
density, the above jump condition shows that dwi/dR is 
almost continuous at the surface [13]. Therefore, in the 
numerical calculation below, we adopt the differentiabil- 
ity condition [23 | 



dwi 
'dR 



0. 



(84) 



At 0(a''^x'^), the evolution equations arc first-order. We 
thus impose the continuity conditions 



[ho{n,)] = 0= [mo (72*)], 
[/J2(72*)] = 0= [fc2(72*)]. 



(85) 



B. Asymptotic Behavior at the Stellar Center 

Before we numerically integrate the interior equations, 
we must understand the asymptotic behavior of the pres- 
sure, density, metric perturbations, and the scalar field at 
the NS center. First, let us focus on the pure GR case. 
In the non-spinning sector, we can use Eqs. ((24)) . ((25)) 
and (P7| to asymptotically expand p, p, M and v about 
i? < i\/ via 



p{R)=Pc + P2R^ + 0{R^), 
2tt 



(86) 



P{R) = - y(Pc + Pc)iPc + iPc)R^ + 0{R^) , (87) 

Air AiT 

M{R) = —p.R' + —P2R' + 0(i?6) , (88) 
3 5 

v{R) = v, + —{p, + 3pc)i?' + 0(i?') {R ^ 0+) , 

(89) 

where Pc, Pc = p{Pc) and are the density, pressure and 
ly at the NS center, respectively. p2 can be expressed in 
terms of pc and Pc through Eq. ([57)) . the TOV equation 
and the EoS. Uc is determined by matching the interior 
solution with the exterior one at the surface. To first- 
order in spin, the solution to Eq. ([32]) that is regular at 
the center is 



* The factor dp/dR is not bounded since p is discontinuous at the 
surface. 



ui{R)=Cjc + —{Pc+Pc)CocR^ + 0{R'') (i?^0+), 

(90) 
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where Wc is a constant that is to be determined by the 
matching condition at the surface. 

Let us now focus on the scalar field and the CS cor- 
rections to the metric perturbation. For the former, the 
solution to Eq. ([55)1 that is regular at the center has the 
asymptotic form 



wi{R) 



Stt 



hoiR) 
moiR) 



12871' 



etc 



1. Solution to C(x'°) 

First, we must obtain the GR solutions at zerotli 
order in spin. There are 4 unknown functions that 
need to be determined, v, A [or equivalently M through 
Eq. (HI])], p and p. The 4 equations that we need to 
solve are Eqs. (gi]), given the EoS. Notice 

that Eqs. (P^ . (P7)) and the EoS form a closed system 
for 7\jf, p and p. We solve these equations as an initial 
value problem using a fourth-order Rungc-Kutta method 
from R = R^ to R ~ TZ^ with the initial conditions 
given in Eqs. ([55|) and ((55)) . where i?e is the core ra- 
dius R^/TZf, <C 1. The radius of the NS can be obtained 
by finding the radius TZ^ where p{TZ*) — and the mass 
of the NS is given by M* = AI (7^* ) . We repeat the cal- 
culation for various pc to obtain a mass-radius relation. 

The remaining equation, Eq. (|25p. can be solved for 
v, again using a fourth-order Runge-Kutta method from 
R = Rg: to TZ^ and with the initial condition given in 
Eq. ([55]) . However, one must be careful that i^(7?.») sat- 
isfies the boundary condition given in Eq. (j71|) . Suppose 
we obtained the trial solution fc'tr(i?) by using the initial 
a{pc+Pc){pc + 3pc)e"'''=/^Wci?c-R^ condition v{R^) = (47r/3)(pc+3pc)i?? (i-e. = 0). Since 



§i{R) ^ d'^R+^i5p,~3p,)d'^R^ + OiR*) (i?^ 0+). 
15 

(91) 

1?^ is a constant that corresponds to dd/dR at the center. 
From Eq. (j39p . the metric perturbation at first order in 
spin has the asymptotic behavior 



{pc +Pc)wc - 167rae'^"/^p2i9c 



+ OiR^) {R^0+). 



i?2 
(92) 



Again, and Wc are to be determined by matching at 
the surface. For the i = mode corrections at quadratic 
order in spin, Eqs. (|42])-(l44l) yield 



0{R^), 

647ra(pc + Pc)e' 
OiR^), 



R' 



uJce 



-i/c/2 



3tt{pc + 3pc) 
hO(i?^) (i?^0+). 



(93) 

(94) 
R 
(95) 



Similarly, Eqs. (gl]) and (|in])-(|lH]) yield the asymptotic 
behaviours for the £ = 2 mode: 



/l2(i?) -/l2c^'+0(i?'), 

1 



k2{R) 



3 



3/i2c - l28TT^a{pc + Pc)e-''-^/'^^'^Ojc 



+8ttI3'0'^^] R^ + 0{R^), 



m2{R) = - g [ih'^c - I28^^a(p,+pc)e-''^'^d',ujc 



(96) 



(97) 



(98) 



47r(/9c + 3pc 



(99) 



Here, is a constant that needs to be determined by 
matching at the NS surface. 



Eq. is shift-invariant, vtr{R) plus a constant Ci, is 
also a solution. The new solution ft,, (i?) -I- Ci, will satisfy 
the correct boundary condition provided 



„l/tr(K.)+C,. _ 



/(7^*): 



(100) 



(101) 



which then yields [58| 

2. Solution to C'(x') 



Next, we solve Eq. p2p . a second-order differential 
equation, for wi. We take advantage of the fact that 
Eq. ([5^ is linear and homogeneous [131 and solve this 
equation with the initial condition given in Eq. (|90l) . We 
choose a specific value for Coc and solve the equation from 
R = Re to 7?.* using a fourth-order Runge-Kutta method 
to obtain the trial solution QtiiR) and the trial moment of 
inertia In calculated from Eq. (|78p . The solution we seek 
is one that satisfies the boundary condition of Eq. (|77p . 
It is improbable that cjtr will satisfy this boundary condi- 
tion. However, due to the homogeneity of the differential 
equation, the product of wtr and a constant is also a solu- 
tion. Hence, we can construct a new solution via CcjCOtr, 
where Cq is a constant, which will satisfy the boundary 
condition, provided that 



C. Numerical Method 

Let us now explain how we solve the interior equa- 
tions numerically to obtain the corresponding interior 
solutions. 



Ci,Wtr(7?.*) =0*1 



~1W 



which then yields 



Co. 



(102) 



(103) 
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Clearly, (like or pc) is a quantity that must be spec- 
ified a priori and controls how rapidly the NS is rotating. 

The scalar field satisfies Eq. ((36)) . which is not homo- 
geneous or shift-invariant. Since the differential equa- 
tion is linear, we take the following approach First, 
we solve Eq. ([36)) with arbitrary values for ■di{R^) and 
'd[{Rf) such that the solution satisfies regularity at the 
NS center [i.e. with an arbitrary value of -d'l in Eq. (|9T)) ] 
to obtain a particular solution in the interior re- 

gion. Next, we look for a homogeneous solution t9)"*(i?) 
in the interior region by solving Eq. ([36]) with a vanishing 
source term on the right-hand-side with arbitrary ■di{R^) 
and "d'liRe) such that satisfies regularity at the 

NS center. The asymptotic behavior of 'd^^^{R) that is 
regular at the NS center is also given by Eq. ([9T)). With 



and at hand, one can construct a generic 

solution to Eq. ([55)) in the interior region that satisfies 
regularity at the NS center as 



(104) 



where is a constant. The constants in Eq. ([?!)) 
and in Eq. ()104p are determined by the matching 
condition at the NS surface given in Eq. ([79| . One can 
solve Eq. ([79)) for and algebraically in terms of 
and 

The CS correction to the gravitomagnetic sector of the 
metric is controlled by Eq. ([39)) . which is a second order 
differential equation for wi. Again, this equation is not 
shift invariant or scale invariant, and thus, we must solve 
it as we did the scalar field above. The matching condi- 
tions at the NS surface are given in Eqs. ([80)) and ([84l) . 
while the asymptotic behavior of the solution that is reg- 
ular at the NS center is given in Eq. ([M)) . The asymptotic 
behavior of the homogeneous solution is given by setting 
a = in Eq. ([92)) . which is the same as Eq. ([90)) with Uc 
being replaced by Wc- 

We have that the solutions for ??i and wi obtained in 
this way are identical to solving Eqs. ([55)) and ([M)) via a 
purely numerical shooting method. 



3. Solution to ©(x'^) 

For the i — mode at 0(a'^x'^), we need to solve 
Eqs. ([43]) and ([44]) to obtain hoc in Eq. ([42]). Since these 
are two coupled first-order differential equations, we can 
solve them as an initial value problem with the initial 
conditions given in Eqs. ()94p and ([M)). Then, we impose 
the continuity condition at the surface to obtain Mqs in 
Eq. ([54]) . The interior solution for Hq can be obtained 
through Eq. ([42]) . By imposing the continuity of Hq at 
the surface, one can determine /iqc- 

For the i = 2 mode, we need to solve Eqs. ^7} and ([ig]) 
for /i2 and fc2- We take the same steps explained above 
for solving Eqs. ([36)) and ([39| . We first obtain the par- 



solutions h™^ and fc™*. The initial conditions at the NS 
center is given by setting a = and /3 = in Eqs. ([M]) 
and ([W]) . As in Eq. ()104p . one can construct generic so- 
lutions for /i2 and k2 that are regular at the NS center 
via 

h]^HR)^h^\R) + C^h^'{R), (105) 
k'riR)^kl^\R) + C'^k\^'{R), (106) 

where Cq is a constant that is to be determined, to- 
gether with another constant Cq in Eqs. (|A1[) and ()A3p . 
with the matching condition at the NS surface given in 
Eq. One can solve Eq. §5^ for Cq and C^ al- 

gebraically in terms of /i2(7?.*) and k2{TZ^). We have 
checked that the solutions obtained in this fashion are 
identical to numericall y solving Eqs. (|47p and 
Riccati method (HMfB. 



via a 



VII. NUMERICAL SOLUTIONS 

In this section, we show the numeri cally obtained NS 
solutions using 4 different EoSs: APR gl, SLy glliil, 
LS220 and Shen 

Wc begin by showing that the results obtained here 
reproduce previously obtained results, but for a wider 
range of EoSs. The left panel of Fig. [3] shows results 
at zeroth order in spin, i.e. the NS mass-radius relation. 
Observe that the SLy curve agrees with [4^ . Notice also 
that the 4 EoSs we adopt in this paper are consistent 
with the existence of the recently found (1.97 ± O.O4)M0 



millisecond pulsar J1614-2230 [42|. The right panel of 
Fig. [3] presents results at linear order in spin, i.e. the 
moment of inertia as a function of mass. The APR curve 
agrees exactly with the results of [i^l)- The left panel 
of Fig. 13] shows i?i/(^2*a//3) as a function of R/TZ^ for 
a NS of mass = 1.4Mq. Observe that -di has its 
maximum near the NS surface. The behavior of these 
curves is consistent with Fig. 5 in j^^]. The top right 
panel of Fig.|3]shows wi/fi* versus R/TZ^,, with 



(107) 



One can see that frame-dragging is strongest near the NS 
center. These curves are consistent with Fig. 6 of (23 |. 
In the bottom right panel of the same figure, we show 
""^1/(^08^^*), which is also peaked in the inner region. 
This behavior is also consistent with Fig. 6 of (23 |. 

Let us now present new results, valid to quadratic or- 
der in spin. Figures [5] and [5] show the metric perturba- 
tions at quadratic order in spin for the £ = and £ = 2 
modes respectively. The interior behavior of ho, and 
^2 is similar to that of The strange behavior of m2 



ticular solutions /i™* and fc^"' with an arbitrary choice of 
/i" in Eqs. ([TO]) and ([^7]) . Next, we look for homogeneous 



We had difficulty in stably carrying out the shooting method due 
to the nearly scale-invariant structure of Eqs. II47I I and 1 148 I I for 
most of the region in parameter space. 
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FIG. 3. (Color online) The mass-radius relation (left) and the moment-of-inertia [Eq. (|76p ] vs. mass relation (right) in GR for 
NSs with various EoSs: APR (sohd), SLy (dotted), LS220 (dashed), and Shen (dot-dashed). 




3 



R/R, R/R. 



FIG. 4. (Color online) Left: Scalar field profile as a function of NS radius with Af* = IAMq and for various EoSs. The 
linear behaviour of i9 near _R — >■ is smooth and differentiable, after multiplying by Pi (cos 0). Right: Correction to the metric 
rotation u) as & function of radius, at linear order in spin, for various EoSs with M, = IAMq. The top plot shows the 0{a'^) 
(GR) contribution loi/Q., [Eq. (fW)) ]. while the bottom plot shows the 0(a'^) (CS) contribution wi/(gcs» ,) [Eq . p8)) ]. In all 
cases, only the I — 1 multipole moment survives at linear order in spin in the small coupling limit. See Sec. IVIBl for details on 
asymptotia. 



in the interior reflects the strange behavior of dp{R)/dR, 
plotted in Fig.[71 The oscillatory structure in these quan- 
tities disappears for a polytropic EoS. The oscillations in 
the EoS, which th en prop agate to 7712, are due to nuclear 
phase transitions [U lisl. I47ll49l[50| . 

Figure [1] already presented p., [Q/{J'^/M^)]{OA/Q, 
6M{ni^^/n^)^{0A/C) and ((5J/J)(0.1/C) as a function 
of M^/Mq for various EoSs. This behavior can be nicely 
fit by 

Ai — exp (fli + biX + CiX^ + diX^ + e^x^) , (108) 
/_ Q 0.1 (5Af rj?,„3 0.1 6J0A\ 

(109) 

with X = M,,/M(7). The estimated numerical coefficients 



are shown in Table HI together with standard errors. Fig- 
ure [1] also shows these fits as a function of M^,/Mq. 
Notice that the CS correction increases the spin angu- 
lar momentum which is consistent with the result re- 
ported in [2^. Moreover, the CS correction increases 
the quadrupole moment while it decreases the observed 
mass. 



VIII. NEUTRON STAR BINARY EVOLUTION 

Let us now study the evolution of a NS binary with 
component masses m^, radii TZ^,a, spin angular mo- 
menta XAfn\S\, dimcnsionless scalar charges Jia and 
quadrupole moment shifts Qa- We use the subscript 
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EoS 


a 


b 


c 


d 


e Error (%) 




APR 
SLy 
LS220 
Shen 


2.90 (0.022) 
2.82 (0.027) 
5.87 (0.071) 
5.35 (0.049) 


-2.30 (0.067) 
-2.31 (0.087) 
-8.07 (0.25) 
-6.86 (0.16) 


1.54 (0.074) 
1.62 (0.10) 
6.60 (0.31) 
5.39 (0.20) 


-0.546 (0.036) 
-0.581 (0.052) 
-2.63 (0.17) 
-2.07 (0.099) 


0.0748 (0.0062) 0.11 
0.0777 (0.0095) 0.12 
0.400 (0.033) 0.49 
0.307 (0.018) 0.36 



[Q/(JVM.)](0.1/C) 



APR -0.802 (0.034) 

SLy -1.79 (0.055) 

LS220 6.72 (0.21) 

Shen 5.01 (0.13) 



-0.675 (0.10) 
1.51 (0.18) 
■16.8 (0.72) 
13.0 (0.42) 



0.578 (0.11) 
-1.37 (0.20) 
15.6 (0.91) 
11.8 (0.49) 



-0.192 (0.053) 
0.675 (0.10) 
-6.66 (0.49) 
-4.94 (0.25) 



0.0150 (0.0091) 0.16 

-0.138 (0.018) 0.23 

1.06 (0.097) 1.2 

0.778 (0.045) 0.67 



-{5M/MQ){nin.s/n,f {0.1/0 



APR 
SLy 
LS220 
Shen 



-2.42 (0.046) 
-3.48 (0.15) 
4.53 (0.24) 
3.00 (0.16) 



-0.503 (0.14) 
2.18 (0.48) 
-14.0 (0.85) 
-10.1 (0.51) 



0.752 (0.15) 
-2.29 (0.55) 
12.3 (1.1) 

8.60 (0.61) 



-0.266 (0.069) 
1.33 (0.27) 
-4.81 (0.58) 
-3.24 (0.30) 



0.0234 (0.012) 0.18 

-0.297 (0.0487) 0.84 

0.673 (0.11) 2.0 

0.445 (0.055) 1.5 



(5J/J)(0.1/C) 



APR 
SLy 
LS220 
Shen 



-1.43 (0.030) 
-0.242 (0.26) 
5.22 (0.11) 
3.94 (0.097) 



-1.73 (0.088) 
-6.30 (0.80) 
-16.1 (0.37) 
-12.9 (0.32) 



1.92 (0.095) 
7.80 (0.91) 
15.9 (0.46) 
12.1 (0.37) 



-0.815 (0.044) 
-3.97 (0.44) 
-7.27 (0.24) 
-5.18 (0.18) 



0.141 (0.0074) 0.13 

0.766 (0.078) 1.0 

1.30 (0.047) 0.56 

0.869 (0.032) 0.98 



TABLE I. Numerical coefhcients for the fitting formula of ft, Q, SM, and 5J as functions of mass using the functional form in 
Eq. (|108[) . Standard errors on each coefficient are in parentheses. For definitions of dimensionless scalar dipole susceptibility 
fi, quadrupole correction Q, mass shift SM, and angular momentum shift SJ, see Eqs. (|69|) . (|66l) . and (|55|) . The last column 
shows the maximum fractional error between values that are obtained numerically and from the fitting formula. 
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FIG. 5. (Color online) i — (left) and 1 = 2 (right) metric perturbations at second order in spin as functions of radius, 
for various EoSs with M* = 1.4M0: ho/{£.cs^*) (top left) and mo/{^cs^*) (bottom left) as defined in Eqs. (|40aP and (gObJ; 
'i2/(Ccsfi^) (top right) and fc2/(^csfi^) (bottom right) as defined in Eqs. (|40ap and (|40cp . 
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R/R. R/R„ 

FIG. 6. £ — 2 metric perturbation at second order in spin as a function of radius, for various EoSs with M» = IAMq: 

Tn2 / {£,csi^*) [Eq. (|40b|) ] for the interior (left) and the exterior (right). Notice the difference in scales. The features in the 
interior correspond to the locations of nuclear phase transitions (see Fig. ^ . 




FIG. 7. Mass density p (top) and radial derivative dp/dR 
(bottom) in the NS interior, for various EoSs with M, ~ 
IAMq. Recall that R is defined by surfaces of constant p, so 
there is no angular dependence to p. Nuclear phase transitions 
are visible as sudden changes in slope. These give rise to the 
features in the metric functions seen in Fig. [B] 



A = (1,2) to denote the A^^ binary component. We also 
use the subscript "12" to denote relative differences. The 
binary orbit can be described by 5 intrinsic and extrinsic 
parameters: semimajor axis a, eccentricity e, the angle of 
periastron w, inclination t, and the angle of the ascend- 
ing node fl (not to be confused with the gravitomagnetic 
metric perturbation w or the NS angular velocity f2*). 

A binary's orbital evolution is affected by dissipative 
effects and conservative effects. The former can be ob- 
tained by calculating the energy flux and the angular 
momentum flux radiated via gravitational and scalar ra- 
diation, which we explain in detail in Appendix |B] In 
this section, we focus on the latter. Here, we are only 



interested in CS corrections on the conservative effect 
relative to the leading GR contribution, and hence, we 
consider CS corrections to Newtonian dynamics. We fol- 
low Gauss' perturbation method [l,[6l,|63] with the con- 
ventions and notation of 0- One must then calculate 
the perturbing accelerations Sa\ and from them the rel- 
ative perturbing acceleration Sa\2 (in any convenient co- 
ordinate system). This vector is then decomposed by 
projecting onto a time-varying orthonormal triad with 
e\ = n\2 and 62 = L"^ (and 63 = 62 x ei so as to complete 
the orthonormal triad). In this triad, the components of 
Sa\2 are defined as 



(110a) 
(110b) 
(110c) 



where inner products are taken with a flat Euclidean met- 
ric. With this decomposition, the osculating orbital ele- 
ments evolve secularly as 









(5a^2e2,i 







he 



■ cost 



(p r)5^ 
he 



sm ( 



f2 cos I 



^ . , y fap 

aS/ft sm H r 

e \ r 



a = 


2a2 
h 




.^e sin 




h 


cos(w + 


0), 


n = 


Wr 

IT 


sin(a; + 


(j)) CSC L 



Here, 



p = a{l~e^) 



(Ula) 
(111b) 

(111c) 

(Hid) 
(llle) 

(112) 



is the semi-latus rectum, r and (f> are the quantities re- 
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latcd to the instantaneous orbital elements, given by 

_ P 



1 



r^-- = h : 



dt 



(113) 
(114) 



with m = mi + TO2 and h is the orbital angular momen- 
tum per unit mass. In all of Eqs. (jllip . the right-hand 
sides are to be orbit averaged as defined in Eq. (|B2p (this 
is appropriate when the time derivatives of the osculating 
elements are much smaller than the orbital timescale and 
there are no resonances). Of course, all of the </> depen- 
dence in #',J5^ and r must be included in the orbit 
averaging. 

In Sec. IVIII Al we calculate the secular evolution from 
the near-zone metric deformation due to the CS correc- 
tion to the quadrupole. In order to calculate the secular 
evolution due to the scalar field, we derive the effective 
dipole interaction in Sec. IVIII Bl With the acceleration 
due to the scalar interaction, we calculate the scalar's 
correction to the binary's secular evolution. 



where ri2 is the separation of the binary components and 
is the unit vector from body 2 to body 1. The ac- 
celeration on body 2 due to body 1 is simply the above 
expression with 1 2. This gives the relative accelera- 
tion 



<ijk> 



15!^i2_ (Q.SiS'l + Q2SiSl) . (116) 



The orbit averages are performed on eccentric orbits 
lying in the x — y plane with the major axis along the 
i-axis, parametrized as [1^ 



x\ = di (cos 4>, sin 0, 0) , 
xl — ~d2 (cos (/), sin 0) , 



di 



7712 



mi 



a{l 



1 + e cos ( 



(117) 
(118) 



(119) 



A. Metric Quadrupole Correction 

Since both the mass monopole and current dipole of 
each body is reabsorbed into physically measured quanti- 
ties, the leading order CS correction is the quadrupole de- 
formation to the metric given by Eq. (j66|) . From Eq. (|65p . 
we can find the acceleration on body 1 due to body 2, 



-15Q2 



<i]k> 
'■12 



oj ok 



(115) 
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[ma(l — e^)] 
d2 



1/2 



(120) 



from which v and all other derived quantities can be cal- 
culated, e.g. n\2 = {x\ —X2)/d= (cos(/), sin0, 0). 

Taking the above expression, decomposing as in 
Eqs. (jllOp . and performing orbit averaging with the pe- 
riod P = 27ra3/2/Mi/2, we find 



n) = 





/^{l - 


e2)2 




3 






/m(l — 


g2)2 




3 






/7^(1 - 


g2)2 



1+2 y'^i,x + ■5*1, yj ^ Si^z cot L ySi^x sinw + Sl^y COSUJ 
QiSi,z {Si_x COS a; — Si_y sinw^ + (1 o 2) , 
QiSi^z CSC L {^Si^x sinw + Si^y coscjj + (1 <-> 2) , 



+ (lo2), (121a) 
(121b) 
(121c) 



and {e)h = = (a)/i, and where the i = x^y,z components of SA,i are taken in the (non-inertial) coordinate system 
where the binary's orbit is in the x — y plane, with pericenter along the -1-5: direction. 

I 



B. Dipole Interaction and Scalar Force Correction density as 



To derive the effective dipole interaction, we start by 
finding an effective interaction Lagrangian between one 
of the body's scalar dipole moments and the scalar field. 
This comes from the cross-interaction part of the kinetic 
term of the scalar field in the action, i.e., by decompos- 
ing the kinetic term of the scalar field in the Lagrangian 



where 



^(a^7?)(9'^79) = /:kin,l + Cvn..2 + Ci^t , (122) 



.-kin, A 



/:kin,int = -/3(9^Z?l)(9^t?2), 



(123) 
(124) 



and £kin,int corresponds to the dipole interaction La- 
grangian density. By substituting t?^ = MA'^^.j/'^A ~ 
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-fi\di{l/rA) into Eq. (|124l) . wc obtain 

ri J \r2 



-Ckin.int = -I3fj.{fl2dij ( I- ) C'ifc ( 



'Pf4f^2di;' { ^) ) , (125) 



where d^^^ denotes the partial derivative with respect to 
x\. 

The dipole interaction Lagrangian can be obtained by 
performing the volume integral of £kin,int as 

= -fiiAiAd^^ j j;^/^- (126) 

By applying Hadamard regularization, as explained 
e.g. in Appendix B of [53| . the above integration can be 
performed to yield 

''12 

= 47r/3— [3(/ii • ni2){n2 ■ n-12) - (mi ' M2)] • (127) 
'^12 

I 



As there are no derivatives on particle locations, this 
gives an effective pairwise interaction potential 



t/int = -L 



int ■ 



This expression agrees with Eq. (5) of [2^. We can find 
the acceleration with (minus) the particle derivative of 
the effective pairwise interaction potential as 



(128) 



From the above we can compute the relative dipole- 
dipole force 



m _ 75 ^cs /^3^3- - 

«12,i - Y^^XlX2C.iG2MlM2 



? 12 / 1^ 



"12 • S2 



1 



(5i • ^2) - 5(7112 • ^i)(ni2 • ^2) 2) 

(129) 



Decomposing aj^2 j inserting this into the Gauss equa- 
tions, and averaging, we find 



256 m4 ( 



l^^~^^)2^i^2MiM2 (^~^ I 2 ('^i^^'^^.x + 5'i,j,S'2,y^ — 5'i,zS'2,2 



cot L 



Si,z ^S'2,2; sincj + 5*2,^, coscj ] 



[n 

and (e)^ = = (d) ,? 



75 Ccs X1X2 



256 m* (1 - e2)' 



,C,3C2^/2i/22 (^) 



7/2 



)]} + (1^2), 
'S'i,z (^'S'2,2; cosw — 52.1, sinoj^ + (1 2) 



75 ^cs X1X2 ^3^3,-,.r,. /"^V^^ 



/ 772 

4n 2^2'^^'^2Ml^2 ~l csct 
■ij 256 TO* (1 — e^)^ V a / 



Si 



,z (52,2; si 



sin a; + S2,y cosoj 



+ (1^2), 



(130a) 
(130b) 

(130c) 



IX. APPLICATIONS TO NS OBSERVATIONS 



A. Massive Millisecond Pulsar J1614-2230 



In this section, we apply the results derived in the pre- 
vious sections to observed NS systems. We consider the 
recently found massive millisecond pulsar J1614-2230 [42| 
and the double binary pulsar system PSR J0737-3039 [3- 
0]. In this section, we consider the maximum CS correc- 
tions allowed within the weak coupling approximation, 
i.e. C = 1- Clearly, such a value for ( violates the small 
coupling approximation. But as we shall see, even with 
such a large ( value, which leads to a strong GR devia- 
tion, NS observations will still not be accurate enough to 
allow for a bound on the theory. 



From measurements of the Shapiro time delay, the 
mass of the pulsar J1614-2230 has been determined to 
be (1.97 ± G.04)Mq [42|. As mentioned previously, the 
effect of a GS modification is to reduce the magnitude 
of the observed mass relative to the GR expectation. In 
Fig. [51 we plotted the mass-radius relation in both GR 
and dynamical GS gravity. We have set the spin pe- 
riod of the NS to be the one observed for J1614-2230, 
P,pi^ = 3.1508076534271(6)ms. When making these 
plots, we did not include 0{a'°x'^) contribution, but we 
have checked that such contributions only affect the re- 
sults to 0.2% at most. As expected, the maximum mass 
decreases as we increase the GS dimensionless coupling 
constant C. Therefore, we can try to place bounds on the 
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Timing parameter 




Pulsar A 


Pulsar B 


Spin frequency 


f. (Hz) 


44.054069392744(2) 


0.36056035506(1) 


Orbital period 


Pb (day) 


0.10225156248(5) 


— 


Eccentricity 


e 


0.0877775(9) 


— 


Projected semimajor axis 


(a/c) sint (s) 


1.415032(1) 


1.5161(16) 


Longitude of periastron 


con 


87.0331(8) 


87.0331 + 180.0 


Advance rate of periastron 


(^) (° /yr) 


16.89947(68) 


[16.96(5)1^ 


Orbital decay rate 


Pb 


-1.252(17) X 10"^^ 




Gravitational redshift 


7 (ms) 


0.3856(26) 




Shapiro delay range 


r (fis) 


6.21(33) 




Shapiro delay shape 


s 


0.99974(-39,+16) 




Inclination 


^(°) 


88.69(-76,+50) 


Mass function 


(Me) 


0.29096571(87) 


0.3579(11) 


Mass ratio 




1.0714(11) 


Mass (Mq) 




1.3381(7) 


1.2489(7) 



" An independent parameter fit of {ui) for pulsar B is consistent with the more precise result of pulsar A Q. 

TABLE IL Timing parameters of the double binary pulsar PSR J0737-3039 0- Measurement uncertainties on last digits are 
shown in parentheses. 



EoS 




APR 


-3.7 X 10"" 


SLy 


-5.1 X 10"* 


LS 


-3.2 X 10"^ 


Shen 


-3.0 X 10"* 



TABLE in. {w)h/{dj)GR for PSR J0737-3039 with various NS 
EoSs, where we set C = 1 (this ratio is linearly proportional 
to C). 



theory by requiring that the maximum observed mass be 
greater than l.QSA/©. Of course, such a method to test 
modified gravity theories is not new; it has been consid- 
ered in e.g. [66{ for Einstein-Aether theory and [63] for 
Einstcin-Dilaton-Gauss-Bonnet theory. 

The horizontal dashed line in Fig. [2] is the lower bound 
on the mass of PSR J1614-2230. For the LS220 EoS, the 
CS maximum NS mass with ^ — 1 becomes less than 
1.93il/Q. This means that we can place a meaningful 
constraint on the theory assuming that this is the correct 
EoS. However, the CS maximum NS mass with C = 1 
with the APR, SLy and Shen EoSs are all above l.gSM©. 
Therefore, whether a meaningful constraint of dynamical 
CS gravity can be placed depends strongly on the EoS. 
Since the EoS has not been observationally measured, 
this degeneracy prevents us from constraining dynamical 
CS gravity with mass-radius relations. 



B. Double Pulsar Binary PSR J0737-3039 

The observed timing parameters of PSR J0737- 
3039 [H-Ql are summarized in Table |lll One can test GR 



from measurements of the post-Keplerian (PK) parame- 
ters [5l| : the (orbital averaged) advance rate of the pe- 
riastron (tb) , the orbital decay rate P, gravitational red- 
shift parameter 7, and the range r and the shape s of the 
Shapiro time delay. 

Let us first look at the CS effect on (w) . Since the spin 
of the secondary pulsar is 100 times smaller than the pri- 
mary one, we only consider CS corrections that depend 
on xi, meaning we only keep {uj)h and neglect {uj)§. By 
taking the ratio of the former to the GR expression for 
(cj)c.R = 3(27r/P)5/3m2/3(l - e^)"! 0,||, we obtain 



Qi 



mr 1 



-Si,z cot L\Si X sin Lo + Si,y cos a; 



(131) 



Since this is proportional to m/a, the CS correction to 
(cO) is of IPN order relative to GR. 

P is proportional to E which is given in Eqs. (jB8[) 
and (|B21[) and both are proportional to {m/ay. By 
comparing this with the GR expectation, which is pro- 
portional to {m/af', the CS correction to P is of 2PN 
relative order. Similarly, one can show that the CS cor- 
rection to 7, r, and s are also of relative 2PN order. 
Therefore, the dominant CS correction to the binary evo- 
lution is in (w) and we can neglect CS corrections to all 
other PK parameters. We can also neglect the CS cor- 
rections to the spin-precession equations, which appear 
at 2PN relative order (see Appendix [C]) . 

In order to place constraints on dynamical CS grav- 
ity, one needs to know the orientation of the spin of the 
primary pulsar since Eq. (jl31l) depends on S\. Unfor- 
tunately, this quantity is currently unconstrained. In 
2004, Ref. [13 attempted to measure this quantity by 
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modeling the intensity variation of pulsar B caused by 
the emission from pulsar A [g^. Doing so, one obtained 
6*5^ ^ 167° ± 10° and 0s, = 246° ± 5°, where 6*5, is the 
angle between SI and the unit orbital angular momen- 
tum vector L', and ^s, is the angle between SI projected 
onto the orbital plane and the direction of the ascending 
node. However, since then, Ref. [i^ found no observable 
major profile changes in the emission of pulsar A, and 
hence the model of Ref. [i^ is ruled out. Nevertheless, 
following Ref. [t^, we adopt the above values for SI to 
give an order of magnitude estimate of the possible mag- 
nitude of the constraint, if this quantity were measured 
in the future. 

By using Kepler's Law = 7Ti(P/27r)^ and the CS 
quadrupole moment deformation found in Sec. I VIII we 
can evaluate Eq. (|13ip for NSs with various EoSs. The 
results for C = 1 ^^rc summarized in Table IIIIl One sees 
that in the small coupling approximation, the ratio of 
\{uj)h/ {uj)gr\ is of C(10~^) at most. For other choice 
of the unit vector SI, |(w);i/(u;)gr| increases by at most 
0.075%. This means that in order to place meaningful 
constraints on dynamical CS gravity from binary pul- 
sar observations, we need to observe (w) and at least 2 
other PK parameters (or the mass functions or the mass 
ratio) within 10^® accuracy (the latter are needed to de- 
termine the masses of the pulsars). With the current 
observational data, we can determine the masses of the 
pulsars using the mass ratio and s, and test dynamical 
CS gravity with (w). However, since the former only have 
accuracies ofO(10~'^), we cannot place any constraint on 
this theory from the double binary pulsar system. This 
conclusion also holds for other binary pulsar systems. 

Notice that |(w)/i/(w)gr| as shown in Table HIT] is 
smaller than |(w)so/('i')GR| ~ 10~^ where (w)so is the 
advance rate of the pcriastron due to GR spin-orbit cou- 
pling and it is given by [t^, l7l| 



(w)gr 



1 4toi -I- 3to2 /i/i/ 



m\ 1/2 
a 



6\/l - e2 mi 
X (2 cos 6s\ + cot i sin 6*5, sin (psi ) , (132) 



where f^-^ is the spin frequency of pulsar A. Since this 
ratio depends on the moment of inertia /, which encodes 
information of the internal structure of the NS, it would 
be difficult to test dynamical CS gravity even if the mea- 
surement accuracies of the PK parameters improved, un- 
less / were determined to high accuracy. 

The results obtained here suggests that pulsar binaries 
cannot currently be used to constrain dynamical CS grav- 
ity. If so, GW observations are the only way to test dy- 
namical CS gravity in the dynamical, strong-field regime, 
as demonstrated in 12911. 



X. CONCLUSIONS AND DISCUSSIONS 

In this paper, we extended the previous work in [23.[58j 
to construct slowly-rotating NSs in dynamical CS gravity 



in the small-coupling and slowly-rotating approximations 
to quadratic order in spin. At this order, we found a neg- 
ative CS mass correction and a positive CS quadrupole 
moment deformation. We applied the former correction 
to test the theory by requiring that the observed maxi- 
mum mass be greater than the lower bound on the mass 
of the massive millisecond pulsar J1614-2230. Unfortu- 
nately, we could not obtain any meaningful constraint 
due to degeneracies with the unknown EoS of nuclear 
matter. Next, we used the quadrupole deformation to 
derive the correction to the evolution of the NS binary. 
Among all PK parameters, we found that the dominant 
CS correction appears in the advance rate of the perias- 
tron (uj). We applied our results to the double binary 
pulsar PSR J0737-3039 0-0 but found that the CS cor- 
rection is too small to be constrained, i.e. it is of relative 
IPN order. This correction is even smaller than the one 
arising from GR spin-orbit coupling, which in turn de- 
pends on the moment of inertia /. Hence, we would not 
be able to constrain the theory unless we knew / to high 
precision. The results obtained here indicate that GW 
observations are the only way to test the theory in the 
strong and dynamical field regime (29j . 

One possible way to extend our work is to construct 
gravitational waveforms for NS binaries. This can be 
done by following the analysis of Ref. [2^. As in the 
BH binary case discussed in the reference, the correc- 
tions to GR waveforms will enter at 2PN order. However, 
one needs to be careful about constructing gravitational 
waveforms for spinning NS binaries in GR, since the spin 
of the NS induces a quadrupole deformation. This defor- 
mation contains information about the internal structure 
of the NS, which also enters at 2PN order [55|- In order 
to obtain 2PN waveforms, one must thus first construct 
slowly-rotating NS solutions in GR to quadratic-order 
in spin. This spin-induced 2PN effect will be strongly 
degenerate with the 2PN CS correction. However, this 
degeneracy might be broken as follows. The tidally in- 
duced quadrupole moment enters at 5PN order (t^ , but 
is enhanced by a factor of 0(7?.* /A/»)^. The tidal effect 
may be detected using future ground-based GW inter- 
ferometers [72| - [76| . Since the internal structure can be 
determined to some extent from tides, the degeneracy 
between the spin-induced term and the CS correction 
might be broken. There is also a GR spin-spin interac- 
tion term at 2PN order in the phase of the gravitational 
waveform. Again, the degeneracies between spins and 
the CS correction term might be broken if the binaries 
are precessing. However, we expect that BH binaries are 
more useful in constraining the theory than NS binaries 
because (i) the degeneracy with the spin-induced term 
would not be broken completely, (ii) the CS deforma- 
tion to the quadrupole moment for a NS is smaller than 
that of a stellar-mass BH due to cancellation between the 
first and the second terms in Eq. ([66)) . and (iii) the spin 
of the NS binaries just before coalescence is expected to 
be rather small [TJ, [zl, [t^, making the CS correction 
even smaller. Also, the radius of a NS is larger than that 
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of a stellar-mass BH, and hence, the latter allows tests 
at smaller length scales, which means that BH binaries 
have a stronger potential to constrain the theory. 

Another possible avenue of future work is to study os- 
cillations of NSs in dynamical CS gravity. By considering 
perturbations of the gravitational, scalar and the matter 
fields, one can investigate e.g. f-, p-, r- and w-modes [78l |. 
The first two modes can be studied under the Cowling ap- 
proximation, where one neglects the perturbations of the 
gravitational and the scalar fields (f- and p-modes in al- 
ternative theories of gravity have been studied in e.g. (tQ^ 
for scalar-tensor theory and [13] for tensor-vector-scalar 
(TeVeS) theory) . This amounts to solving the perturbed 
matter equation V^STJ^J"^ = 0, where STJ^J"^ is the per- 
turbed matter stress-energy tensor. At zeroth order in 
spin, since the background is exactly the same as GR, we 
expect that f- and p-modes to be identical to those in GR. 
To study w-modcs, one needs to include gravitational 
and scalar field perturbations. W-modes in alternative 
theories of gravity have been investigated in e.g. |8l| for 



scalar-tensor theory and [82[ for TeVeS. In these theories, 
one finds that axial metric perturbations decouple from 
scalar and polar metric ones. Hence, w-modes can be 
studied by solving the axial perturbation equation as an 
eigenvalue problem. In dynamical CS gravity, we expect 
that polar perturbations will decouple from the rest, as in 
the case for BHs [83l - [85l |. However, again, since the back- 
ground is the same as in GR, the polar perturbation equa- 
tions should be identical to those in GR. At linear order 
in spin, one can look at r-modes by studying the toroidal 
oscillation of the NS in the Cowling approximation. As in 
GR [1^, 113] and in TcVeS [s^ , the spectrum for toroidal 
oscillations of slowly-rotating stars should also be con- 



tinuous in dynamical CS gravity. However, this might 
be an artifact of the slow-rotation limit, since in GR, the 
spectrum becomes discrete for fast- rotating NSs [1^ [qO] ■ 
In principle, if the maximum and/or the minimum of the 
toroidal oscillation frequencies and the moment of inertia 
are independently determined, one can break the degen- 
eracy between the gravitational theory and the EoS, and 
hence test GR [S^ . This looks challenging from an obser- 
vational point of view. At quadratic order in spin, there 
can be corrections to f-, p-, and w-modcs, but obviously, 
they are suppressed by 0{C,y^). 
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Appendix A: The £ = 2 Mode Exterior Solutions at Second Order in Spin 
We solve Eqs. (j46|)- (|48| and obtain the t = 2 mode exterior solution at 0(a'^x'^): 
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where Cq is another integration constant. One might 
think that the above metric components diverge at in- 
finity. However, one can show that hf^^ = 0{R^^), 
mf^ = 0{R-^) and kf^ = 0{R-^) as /? oo. 



Appendix B: Energy Flux 
and Angular Momentum Flux 

In this appendix, we look at the dissipative effects of 
the binary evolution, i.e. the energy and angular mo- 
mentum fluxes radiated via gravitational and scalar ra- 
diation. 



1. Energy Flux 

The rate of change of the orbital binding energy Ei, 
must be balanced by the energy flux J- carried away 



from the system by all propagating degrees of freedom. 
In dynamical CS gravity, there are two such quantities, 
the metric perturbation (h) and the scalar field and 

thus Eb = -T(h) - ^{-d) = E^''^ + E^'^'^- Each of these 
contributions can be split into a GR term plus a CS 
term, noting of course that E^J = 0, we then have 
Eb = eW + 5eO^) + dE^'^\ where the S's are to remind 
us that these are CS corrections. 

In either case, for any field (p with stress-energy tensor 
T(^), the energy flux is (see Sec. VI of (H) 

= lim / iT^fn') r^dQ, (Bl) 

r-i-oo Jg2 \ / LU 

where the orbit average of any quantity Q is defined as 

where T is the orbital period, </> is the orbital phase, and 
the Jacobian (f> must of course be included. 
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a. Scalar Field 

Let us consider the energy flux associated with the CS 
scalar field The scalar stress-energy tensor was given 
in Eq. ([S])- The radiative far-zone solution for the scalar 
field ^''^ is M, 



(B3) 



where iiij is the magnetic-type quadrupolc tensor of the 
source, 



X2 fJ-2 



(B4) 



We could rewrite these expressions in terms of the struc- 
ture constants of Eq. (|B7|) . but the resulting expression is 
rather long and unilluminating. It may seem surprising 
at first that the different components of A* play unequal 
roles, but the binary orbit sets up a preferred coordinate 
system that treats components of A* differently. 

b. Metric Perturbation 



The leading-order correction to the effective GW 
stress-energy tensor is given by (see and Sec. VI 
of 11) 



There are other contributions to -d^^, but they are sup- 
pressed by the ratio of the orbital and radiation-reaction 
or precession tiniescales, so we neglect them here. 

Inserting this far-zone solution into T^'^'> and this into 
the energy flux formula, we obtain 



47r 
15 



13 



(B5) 



LTT t.lj 



(B12) 



where ()a is a short- wavelength average, hij is the metric 
perturbation in GR, is the correction to the metric 
perturbation due to CS gravity, and TT stands for the 
transverse-traceless projection. 



which upon expansion returns 



768 m4 



{A\,l2 + 2(A • v,2f 



+ 3A2(ni2 • vi2 f - 12(ni2 • wi2)(A • ni2)(A • V12) 
+ 18(A-ni2f(ni2-«i2)'))^ . (B6) 

We have here defined 



(B7) 



A 



ij^kl 



P^kPJl 



-P^jPku (B13) 



with Pij ~ Sij 



j the projector onto the plane perpen- 
dicular to the line from the source to a FZ field point. 
The leading-order expressions for hij and in the far- 
zone are (see Eq. (118) of [H) 



so as to agree with the definition of A* in [53[ when the 
two bodies are BHs. In the circular orbit limit, (nj2 ' 
V12) and our Eq. ((B6)) agrees with Eq. fl37) of [ssj. 

By performing the orbit average, we find that the cor- 
rection to the rate of change to the orbital binding energy 
due to the scalar field is 



JijW ^ _ 2A2/i(e) + 2Ai/2(e) + A^/gj, 



where 
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(B15) 



+ 3ni^2 [5 {n\2^ilk) (ni^mi) - MifcA4] } + (1 ^ 2) , 



(B8) where lij = mixjx] + (1 2) is the mass quadrupole 



tensor. 



The effective stress-energy tensor of Eq. (jB12|) must be 
inserted into Eq. (|Bip to calculate the correction to the 
rate of change of the orbital binding energy due to the 
GS correction to the metric perturbation. Inserting the 
previous expressions and performing the integral over dfl 
gives 
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Similar to the scalar field case, in the circular orbit, 
Eq. (|bT6| divided by the GR energy flux -Egr agrees with 
Eq. (142) of [11]. Performing the orbit-averaging yields 
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(B20) 



Rewriting the energy flux in terms of the dimensionless 
structure constants, we find 
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before [Eq. (|B3p ] into the stress energy tensor [Eq. (O], 
and inserting this into Eq. (jB22p . we have 
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Expanding in terms of time derivatives of the quadrupole 
tensor gives 
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Finally, performing the orbit-averaging as before, with z 
lying perpendicular to the orbital plane, gives 
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2. Angular momentum flux 

In this appendix, we present the rate of change of the 
orbital angular momentum, as induced by the propagat- 
ing scalar field and metric perturbation. 



a. Scalar Field 

The angular momentum flux for the scalar field 
which has stress-energy tensor T*^'') is [6^ 



The angular momentum flux associated with the met- 
ric perturbation at null infinity can be written as (see 
Eq. (4.22') in (Hi) 



7 ttJ TT 2 
"pa aq 



-n h^^ h'--r 
2^P'^ab,q'T'ab ' 



TT ;.TT„3\ 

I X.UJ 

(B27) 



This corrects a well-known mistake in the paper of Pe- 
ters [13 . If the metric perturbation consists of a GR 
term plus a CS correction, we can then calculate the CS 
modification to the angular momentum fiux via 



SL 



-eijk lim 

] — >oc 

-eijk lim 



Si 



Si 



r^dn 



(B22) 
(B23) 



Only the parts of T^; which decay as r""^ contribute a 
finite part. One can verify that the parts that decay as 
actually vanish identically prior to taking the limit 
to spatial infinity. Inserting the same far-zone solution as 



6E 



iDTT r^oo 



TT 

■9 



{y^b,qKib + ^lib,q^Zb) 



I^TT r-TT 

"pa "aq 

dfl . (B28) 



As before, only the parts of the integrand that decay 
as and in the first and second terms respectively 
contribute a finite part. One can verify that any seem- 
ingly divergent terms actually vanish prior to taking the 



26 



limit to spatial infinity. Inserting the same far-zone solu- 
tion as before [Eqs. (|B14|) - (|B15p ] into the above expres- 
sions, inserting the expressions for the time dependence 
of ^ij and orbit averaging, we finally obtain 

(B29) 



16 TO^ 

where we have defined 



1 + f ^2 + le/) (^S^S^ + 5f 5|) , (B30) 



9 2 1 

r +4 
1 



L={1+ "-e" + 4eM ( SlSl + S\S'^ ] , (B31) 



1 - ) S'iS^ + ( 1 + 6e^ -I- -e* ) 5f S| 



2 ( 1 + 3e2 + ^4 ) SIS^ 



(B32) 



One can, of course, check that in the zero eccentricity 
limit 5E^ = ^5L\, where Vl = 2t:/T is the orbital fre- 
quency. 



Appendix C: Spin Precession due to 
Scalar Interaction 

Using the relationship in Eq. ([70|) . we can rewrite the 
scalar dipole-dipole interaction as an additional spin-spin 
interaction. 



6C 



ss 



25 ^cs MlM2^3^3 



256 TO^ 



[3(5i-ni2)(52-ni2)-(5i-52)] , (CI) 



is of the same form with only a different prefactor. This 
interaction will lead to additional spin precession (to be 
added to that already present from GR) of the form 



5S[ 



ss 



25 ^cs M1M2 ^3^3 



256 m4 



3(ni2 • S2)n- 



(C2) 



and similarly for body 2 with 1 <-> 2. 

The dipole-dipole interaction does not modify the spin- 
orbit coupling at all, so the orbital angular momentum 

does not appear in the above expression as it does in 
GR. We expect that the 0{a'^) conservative correction 
to L' appears at 2PN order or higher, so the correction 
to the orbit-induced precession will only appear at higher 
than 2 PN order. 

This effect is interesting in that it comes in at the same 
PN order as in GR. However, currently, there are no bi- 
nary systems with sufficiently well-modeled spin preces- 
sion that could be used to measure or place constraints 
on Chern-Simons gravity. Indeed, the spin precession 
periods for currently known pulsar binary systems arc of 
the order of hundreds of years [11]. Therefore, it seems 
unlikely that this phenomenon will soon be used for con- 
straining modified theories. 

Precession is also caused by a monopole-quadrupole 
interaction [9^, llOCf . which enters at the same PN or- 
der as the spin-spin interaction [HH]- Since there is a 
quadrupole moment deformation in dynamical CS grav- 
ity, there would be a CS correction to the monopole- 
quadrupole interaction of the form 



{S- 



' 12 



* — 2 — 3 — 



ri2 



L-SAe,,kUSt (C3) 



where Ca is the compactness of body A. Compare this 
with the leading spin-spin interaction present in GR, 
e.g. Eq. (5b) of [93] • It is clear that this contribution 



and similarly for body 2. Notice that this is of the same 
PN order as the correction to the spin-spin interaction 
shown in Eq. (|C2]) . 
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